import pydicom
import os
import numpy as np
import datetime

imagedir = r"C:\Users\jerry\Documents\WXWork\1688850228358093\Cache\File\2020-04\PAT0005\S0001"
datasets = tuple([pydicom.dcmread(os.path.join(imagedir, filename)) for filename in os.listdir(imagedir)])
# skip files with no SliceLocation (eg scout views)
slices = []
skipcount = 0

for f in datasets:
    if hasattr(f, 'SliceLocation'):
        slices.append(f)
    else:
        skipcount = skipcount + 1

print("skipped, no SliceLocation: {}".format(skipcount))
# ensure they are in the correct order
slices = sorted(slices, key=lambda s: s.SliceLocation)

# create 3D array
img_shape = list(slices[0].pixel_array.shape)
img_shape.append(len(slices))
img3d = np.zeros(img_shape)

# fill 3D array with the images from the files
for i, s in enumerate(slices):
    img2d = s.pixel_array
    img3d[:, :, i] = (img2d * s.RescaleSlope + s.RescaleIntercept)

maxvalue = np.max(img3d)
print(maxvalue)
print(np.where(img3d == maxvalue))

def gettimeformat(timevalue):
    return '%H%M%S.%f' if timevalue.find('.') > 0 else '%H%M%S'

def getdatetime(datestr, timestr):
    date = datetime.datetime.strptime(datestr, '%Y%m%d')
    time = datetime.datetime.strptime(timestr, gettimeformat(timestr))
    return datetime.datetime.combine(date, time.time())

def getserialdatetime(ds):
    return getdatetime(ds[0x0008, 0x0021].value, ds[0x0008, 0x0031].value)

def getacqdatetime(ds):
    return getdatetime(ds[0x0008, 0x0022].value, ds[0x0008, 0x0032].value)

ds0 = slices[0]
halflife = 6586.2
decayTime = getserialdatetime(ds0) - getdatetime("20200410","085100") #getserialdatetime(ds0)
injectedDose = ds0.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose
decayedDose = injectedDose * pow(2,-decayTime.seconds/halflife)
weight = ds0.PatientWeight

SUVbwScaleFactor = (weight * 1000 / decayedDose)
print(maxvalue * SUVbwScaleFactor)